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ABSTRACT 

We have developed a parallel, simple, and fast hydrodynamics code for multi- 
dimensional, self-gravitating, adiabatic flows. Our primary motivation is the study of 
the non-linear evolution of white dwarf oscillations excited via tidal resonances, typ- 
ically over hundreds of stellar dynamical times. Consequently, we require long term 
stability, low diffusivity, and high algorithmic efficiency. An Eulerian finite-difference 
scheme on a regular Cartesian grid fulfills these requirements. It provides uniform 
resolution throughout the flow, as well as simplifying the computation of the self- 
gravitational potential, which is done via spectral methods. In this paper, we describe 
the numerical scheme and present the results of some diagnostic problems. We also 
demonstrate the stability of a cold white dwarf in three dimensions over hundreds of 
dynamical times. Finally, we compare the results of the numerical scheme to the linear 
theory of adiabatic oscillations, finding numerical quality factors on the order of 6000, 
and excellent agreement with the oscillation frequency obtained by the linear analysis. 
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>^ ■ 1 INTRODUCTION 



White dwarfs (WDs) are known to be common endpoints of 
stellar evolution. A significant amount of evidence suggests 
that both stellar mass black holes (< 10^ M0 ) and neutron 
stars are also rela tively common. More recently, both the- 
oretical ( see, e.g. , iMadau fc Reesl l200j) and observ ational 
(see, e.g., IColbert fc PtakI 120021: iGerssen et alJl2002^ stud- 
ies have implied the existence of intermediate mass black 
holes (lO'^-lO^M© ). As a result, it appears inevitable that 
white dwarf-compact object binaries will form. This may be 
especially likely within cluster environments. 

After formation, the subsequent evolution of a white 
dwarf-compact object binary will typically be driven by 
gravitational radiation. As the system passes through res- 
onances between the normal mode frequencies of the white 
dwarf and harmonics of the orbital frequency, it is possible 
to resonantly excite oscillations on the white dwarf. Even 
small amounts of energy transfer may have a non-negligible 
impact upon the orbit, possibly with consequences for grav- 
itational wave detections of such systems (e.g., by LISA). 
Large energy transfers may result in heating and, possibly, 
the detonation of the white dwarf, leading to an exotic type 
I supernova and, perhaps, a subsequent 7-ray burst. In order 
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to assess the magnitude and likelihood of such scenarios, it is 
necessary to understand the mode excitation process in de- 
tail. For the linear regime, this has been done iRathore et alJ 
I2OO3I |20o3), and it was found that, depending upon the 
initial conditions, it is possible to excite modes with large 
enough amplitudes that the validity of the linear theory be- 
comes questionable. Therefore, it is necessary to investigate 
the mode evolution in the non-linear regime. This is most 
directly done via numerical hydrodynamics simulations. 

A number of hydrodynamics codes which may be used 
for this purpose curre ntly exist. Two such codes, ZEUS 
llStone fc Normanll992l) and Flash iFrvxell et al.l2000l) have 
been developed to be generic hydrodynamic engines. Such 
codes provide access to a sophisticated suite of hydrody- 
namic simulation tools. However, they also have the disad- 
vantage of being complicated to use and, perhaps, subop- 
timal for our specific problem. In addition, to a good ap- 
proximation, the white dwarf oscillations are adiabatic, and 
hence detailed treatment of shocks and entropy generation 
are unnecessary. 

iMotl et al.l i2002t) have developed an adiabatic hydro- 
dynamics code, primarily for studying binary mass transfer. 
However, the choice of a cylindrical grid, while useful for 
the mass transfer application, is problematic for the case 
of a pulsating white dwarf, where it is important to main- 
tain uniform resolution throughout the star. Furthermore, a 
cylindrical coordinate system complicates the numerical ad- 
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vection scheme. These difficulties are avoided with a Carte- 
sian grid, an additional advantage of which is that the Pois- 
son equation can be solved easily and efficiently via spectral 
methods. 

We present a simple hydrodynamics code with some di- 
agnostics and an example application. This is done in seven 
sections with ^|21 reviewing the hydrodynamic equations, 
outlining the differencing scheme used, ^ describing the 
method used to solve the Poisson equation, JSl presenting 
some tests of the code, JUa-PPlying the code to an oscillat- 
ing white dwarf, and i|7| containing concluding remarks. 



2 GOVERNING HYDRODYNAMIC 
EQUATIONS 

There is considerable freedom in the choice of macroscopic 
quantities used to describe ffuid flows. Our choice was pri- 
marily dictated by the numerical convenience of the sourced 
advective form of the hydrodynamic equations. In addition, 
since we are restricting ourselves to adiabatic flows, it is 
convenient to use the entropy rather than the energy as a 
thermodynamic variable. We therefore chose the following 
five quantities to describe the fluid flow: mass density (p), 
entropy density (s), and momentum density (J). 

The equations for p and s have a purely advective form. 



ds „ 



vs = 0, 



(1) 
(2) 



which correspond to the conservation of mass and entropy.^ 
The equation for J can be written in a sourced advective 
form, 



-VP - pV^ + { , (3) 
where the pressure (P) is given by an equation of state. 



(4) 



the self-gravitational potential ($) is determined by the 
Poisson equation. 



= 47rGp, 



(5) 



and f is any additional external force per unit volume act- 
ing on the fluid (e.g., an external gravitational field and/or 
Coriolis forces). 



3 DIFFERENCING SCHEME 

In one dimension, the use of a staggered mesh avoids the 
interpolation of the flow velocities to the cell boundaries. 
With a zone-centred grid, the velocities would have to be 
interpolated, which would complicate the advection step in 



^ Note that s is the entropy per unit volume and not tlie specific 
entropy. Hence, in our notation, the adiabatic condition is 

dt \p) 

where d/dt is the convective derivative. 
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Figure 1. The geometry of a zone-centred, uniform Cartesian 
grid is shown. Here, A can be any of the five evolved quantities 
(p, s, and J) or the gravitational potential (#). 



the momentum conservation equation However, in mul- 
tiple dimensions, the interpolation of vector quantities (e.g., 
the momentum density) cannot be avoided by the use of a 
staggered mesh. Therefore, we use the conceptually simpler 
zone-centred grid. 

Casting the hydrodynamic equations in a sourced ad- 
vective form allows the explicit conservation of mass, en- 
tropy, and momentum (insofar as the source terms allow). 



3.1 Advection 

The advection steps in equations QISJl may be integrated 
to yield flnite-difference equations for a given cell 



AA: 



At 
AV 



(6) 



where AA is the change in the quantity A due to fluid ad- 
vection, At is the time step, AV^ is the cell volume, A±i are 
the fluxes of the quantity A at the ±ith boundary of the 
cell, and AS'i is the area of the cell surface normal to the 
ith direction. 

In general some interpolation is required to determine 
the values of the fluxes at the boundaries of the cell. We 
break the interpolation of the fluxes into an interpolation 
over the fluid velocity and an interpolation over the advected 
quantities. 



A±i = A±iU±i 



(7) 



where v±i is the interpolated component of the velocity nor- 
mal to the ±ith cell face at the cell boundary, and A* is the 
interpolated value of the advected quantity. The v±i are de- 
fined by 



(8) 
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where ii* and v±i are the values of the fluid velocity in the 
ith direction at the centre of the current cell, and the centres 
of the neighbouring cells in the ±ith directions, respectively. 

A numerical difficulty with the interpolation of the ad- 
vected quantities is that advecting the volumetric densities 
tends to generate unphysically high velocities in low cells 
with low mass density. We circumvent this prob lem by us- 
ing consistent transport jStone fc No rman* 1992), in which 



it is the specific quantities that are interpolated, i.e. 



(9) 



where and {\/p)±i are the interpolated values of p and 
the specific quantity X/p at the ±ith boundary of the cell. 

The choice of the method used for interpolating the ad- 
vected quantities has to be made carefully, so as to avoid 
introducing instabilities in the finite-difference scheme. Sev- 
eral such methods exist, of which we have chosen to use up- 
winding methods. These methods provide stability by clip- 
ping new local extrema, and limit diffusivity by interpolat- 
ing quantities to the boundary in a way that accounts for 
the difference between the velocities associated with the up- 
wind and downwind characteristics. Upwinding methods of 
various orders exist, with the the higher-order methods be- 
ing necessarily more computationally expensive. The three 
methods we have implemented are the donor cell (zeroth- 
order), van Leer (first-order), and piecewise parabolic ad- 
vection (PPA; second-order) methods. 

3.1.1 Donor Cell Upwinding 

The donor cell method is a zeroth-order upwinding scheme, 
approximating the spatial distribution of a given quantity, 
g, as a step function. In this method, all information from 
the downwind cell is ignored, i.e. at the —ith cell boundary 



l-i = 



ff v^i > 
if v-i < 



(10) 



For a given cell, this only requires information from the near- 
est neighbours. In practise, donor cell upwinding is highly 
diffusive (see, e.g., H5.H . and hence was not used beyond the 
testing stage. 

3.1.2 van Leer Upwinding 

The van Leer upwinding metho d is a fir st-order method 
first described by its namesake (Ivan Leer 1977a b. 1979). 
In contrast to the donor cell method, the distribution of q 
is approximated by a piecewise linear function. The slopes 
of these linear functions are given by the so-called van Leer 
slopes, defined below for a given cell along the ith direction, 



! '2{q+i - q){q - q-^) 
Ax'{q+, - q^i) 




if (<?+»- g)('7 - 5-0 > 



otherwise 

(11) 

In terms of the van Leer slopes, the upwinded value of the 
quantity q at the —ith cell boundary is given by 



q-i + ^ (^Aa;' - v^iA.t^ dq^i if V-i ^ 
g - i (^Ax" + v-iAt^ dq^ if «_i < 



(12) 



where the notation dg±j denotes the van Leer slope in the 
ith direction for the neighbouring cell in the ±jth direction. 
The van Leer method prevents the introduction of new lo- 
cal extrema, and hence ensures stability in the advection 
scheme. When the van Leer slopes vanish, the scheme re- 
duces to the donor cell method. Note that, because van Leer 
upwinding uses the van Leer slopes of neighbouring cells, it 
requires information from both the nearest and next-nearest 
neighbours. 



3.1.3 PPA Upwinding 

The PPA method i s a second-order upwinding m ethod orig- 
inally developed bv lColella fc WoodwardI lll984H . It approxi- 
mates the distribution of g by a piecewise parabolic function. 
The essence of the method is the determination of the mono- 
tonized left and right interface values, qL and gn, which are 
computed via equations (1.6)-(1.10) in.Colclla fc Woodw ardI 
|1984). In terms of qL and qn, the upwinded value of q at 
the —ith cell boundary is given by 



qR,-i + i{q-i - qn-i) 
+ ^{1 ^ Oi2q-^ - qR.-^ ' 

qL +S,{q - qL) 



QL-i) 



if v-i ^ 



if V-i < 



(13) 

where ^ — V-iAt/ Ax^ . This requires information from the 
nearest three neighbours. 

The PPA method is substantially less diffusive than the 
van Leer method. This is especially notable at discontinu- 
ities, where the profiles generated by PPA are significantly 
steeper than those generated by the van Leer scheme. How- 
ever, the improvement comes with a relatively high compu - 
tational cost. It has been found bv lStone fc NormanI |^9^ 
that, typically, increasing the grid resolution is a computa- 
tionally more efficient way to obtain greater accuracy. For 
this reason, unless explicitly stated otherwise, we use the 
van Leer upwinding method. 



3.2 Artificial Viscosity 

In Eulerian upwinding schemes, shocks can lead to numerical 
instabilities. If resolving shocks is critical, the instabilities 
may be cured via the introduction of Riemann solvers (capa- 
ble of localising a shock to a single cell boundary). However, 
if resolving shocks is unnecessary, it is significantly easier to 
introduce an artificial numerical viscosity to smooth them 
out. Several prescriptions for implementing numerical vis- 
cosity can be found in the literature; we chose to implement 
the von Neumann-Richtmyer scheme because of its ability to 
produce the correct shock propagation velocity and its low 
dissipation far from shocks (a direc t result of the fact that i t 
acts only in regions of compression; IStone fc Normanlll992ll . 
This scheme takes the form of defining a viscous pseudo- 
pressure for each direction which is non- vanishing in regions 
of compression only: 



;2 

I p 



dx^ 



if 



< 



dv' 
dx^ 

otherwise 



(14) 
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for i = X, y, z, where I is the length scale over which shocks 
are to be smoothed. The associated source term for equation 
iPJ is given by 



Vise 



(15) 



Typically this will smooth a shock front over a number of 
cells-a distance that is usually much larger than the natural 
shock depth. It should also be noted that a strictly correct 
treatment of shocks is precluded by the adiabatic condition, 
equation ^ . This can be remedied by the inclusion of a vis- 
cous source term in the entropy equation. However, since for 
the applications we envision shocks will result in the rapid 
thermalisation of the kinetic energy of the stellar oscilla- 
tions, their mere production may make a purely hydrody- 
namic description inapplicable. In particular, thermonuclear 
processes could dominate at such a point, and thus neither 
the added complexity and computational overhead of the 
Reimann solver methods nor the complication of an entropy 
source term are required. 

3.3 Momentum Source Terms 

In addition to advection, the momentum density evolves 
due to pressure gradients, self-gravity, and external forces 
(if any). We have found that simply finite-differencing VP 
leads to a less stable system than calculating the gradient 
via partial derivatives of the equation of state, and finite- 
differencing in p and s. In contrast, the gravitational ac- 
celeration is obtained directly in terms of a second-order, 
finite-difference of the gravitational potential (the details of 
solving for which are presented in Q . The finite differencing 
of the viscous force is performed in two steps: (i) determin- 
ing the viscous pseudo-pressure, and (ii) finite differencing 
the viscous pseudo-pressure to obtain the viscous force di- 
rectly. In finite difference form, the viscous pseudo-pressure 
is defined by 



P±i+ P f ""'±1 



< 



Q±> = <j 2 V Ax> 

otherwise 

(16) 

for i — X, y, z. The dimensionless coefficient rj is approxi- 
mately the number of cells over which discontinuities are to 
be smoothed. Typically, we find = 2 to be adequate. The 
viscous force is then determined by 



^ VIS' 



Aa;' 



(17) 



Therefore, excluding external forces, the source terms in 
equation are given by 

dP\ p+i-p-i f dP\ s+, -s_i 



dp J 2Ax' \ ds I 2Axi 

for i = X, y, z. 

When using a barotropic equation of state, P{p), it can 
be convenient to write the source terms in terms of the spe- 
cific enthalpy, h, 

'h+t - h-i $+, - 



2Ax^ 



+ 



2Axi 



+ F ■ 

~ VIS' 



(19) 



for i — X, y, z. An example of when this is useful will be 
discussed in Note that in this case, the entropy equation 
is superfluous. 

3.4 Courant-Friedrichs-Lewy Time Step 

The stability of our explicit finite-difference scheme requires 
that the time step should satisfy the Courant-Friedrichs- 
Lewy (CFL) criterion. This corresponds to the physical con- 
sideration that, in a single time step, information should 
only propagate into a given cell from the neighbouring cells 
which are used to compute spatial derivatives at that point. 
A time step that is too large would require information from 
more distant cells, which is not available in the differencing 
scheme. Therefore, for stability. 



At SC tcFL , 

where the CFL time is defined by 
Ax Ay 



iCFL 



Az 



Cs + \V^\ ' Cs + \vy\' Cs + 



(20) 



(21) 



where Cs i s the local adiabatic sound speed (see e.g. 
iMotl et alJ l2002l : IStone fc NormanI Il99l and references 
therein). In addition, the inclusion of an artificial viscos- 
ity imposes the additional requirement that the time step 
does not exceed the timescale for diffusion across cell width 
length-scales: 



/ Ax 



Ay 



Az 



(22) 



\477|Au^| ' 4:ri\Avy \ ' 4r?|A?;^ 

(see e.g. lStone fc Normanll993) . In practise, for many oper- 
ator split methods, taking the time step to be the CFL time 
does not ensure stability. Rather, it is necessary to take At 
to be some fraction of tcFL or tvisc. In practise, we find that 
a robust choice is 



At < -min(tcFL,tvis. 



(23) 



From equation H21|l it is clear that the cells with the 
highest velocities (both kinetic and sound) will provide the 
most stringent limits on the time step. An example is the 
case of cells constituting the vacuum surrounding a star. 
In practise, for numerical reasons, no portion of the grid 
can have vanishing mass density. Therefore, we take 'zero' 
density to be some small fraction (typically, 10~*) of the ini- 
tial maximum density. As a result, the vacuum is physically 
insigniflcant. Nonetheless, because of their large accretion 
velocities (though negligible momentum densities), the vac- 
uum cells can be the limiting factor in determining the time 
step. To avoid this problem, we impose a velocity cap, so 
that the CFL time is set by only considering cells with den- 
sities larger than, say, 10~® of the maximum density.^ The 
remaining cells have their velocities capped at 



Ax Ay 
'At' 'At' 



Az 
'At 



(24) 



so as to not drive the time step down. While this explicitly 
violates the hydrodynamic equations presented in it does 
so in a physically negligible manner. 

^ What is important is that the density cut-off used for the CFL 
time is large enough to exclude the vacuum cells. 
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We use operator splitting to separate the source and ad- 
vection contributions to the evolution of the fluid quantities 
at each time step. However, we do not use directional split- 
ting, making our scheme a variation of the unsplit method of 
van Leer. Thus, a single time step is taken in three stages: 
(1) taking half of the source step, (2) performing the up- 
dates due to advection, and (3) taking the second half of 
the source step. The gravitational potential is calculated at 
each source sub-step. 

3.5 Boundary Conditions 

Because the upwinding methods require information about 
neighbouring cells, it is necessary to provide a boundary 
of ghost cells along the outer edges of the grid. As these 
ghost cells are not evolved themselves, they require some 
prescription for assigning the evolved quantities to them. 
We have implemented three types of boundary conditions: 
fixed, replicated, and outflow. 

The first, and simplest, is the fixed boundary condition. 
In this prescription, the boundary cells are fixed to have 
'zero' density, entropy density, and momentum flux. This 
tends to limit the velocity of the 'zero' density vacuum by 
not providing a boundary momentum flux. 

The second set of boundary conditions consists of repli- 
cating the last set of cells in the grid. This provides a slightly 
more realistic set of boundary conditions, allowing the accre- 
tion of the 'zero' density vacuum to stabilise through hydro- 
dynamic balance. However, if a physically significant portion 
of the fiow is crossing the boundary, then this is significantly 
superior to the first scheme. 

The third set of boundary conditions implemented are 
the so-called outflow boundary conditions. In this prescrip- 
tion, fluid is allowed to flow off the grid but not into it. 
In order to prevent the boundaries from physically affect- 
ing the fluid on the grid, the boundary values for density 
and entropy are chosen to preserve hydrostatic equilibrium 
in the last grid zone. Note that this does not stop the fluid 
from advecting off the grid through this zone. As a result, 
this will minimise the creation of spurious reflections at the 
boundaries. For a self-gravitating fluid configuration that is 
initially contained entirely within the grid, this provides the 
most realistic set of conditions. 

3.6 Parallelisation 

The primary purpose for the development of our code is to 
perform high resolution studies of the non-linear evolution of 
normal modes in white dwarfs. The resulting computational 
requirements necessitate high-performance computing. Be- 
cause the sourced advection step for a given cell depends 
only upon cells in its immediate neighbourhood, it naturally 
lends itself to a straightforward parallelisation scheme. This 
takes the form of dividing the entire grid into a number of 
sub-domains, each of which are handled by a separate pro- 
cess. Because interprocess communication incurs substantial 
performance penalties, we need to choose a domain decom- 
position that minimises the communication required. The 
source of interprocess communication in each sourced advec- 
tion step is the need for neighbour data around the edges of 
each sub-domain. Therefore, the time penalties due to inter- 
process communication are dictated by the surface area of 



each sub-domain, as well as the depth of neighbours that is 
necessary (one for donor cell upwinding, two for van Leer up- 
winding, and three for PPA upwinding). Hence, minimising 
the surface area of each sub-domain minimises the interpro- 
cess communication. 

We have chosen to implement our code in the C-|— I- 
programming language. This choice is motivated by consid- 
erations such as modularity of design, flexibility, efficiency, 
ease of code reuse, and extensibility. For example, using the 
object-oriented paradigm in the C-|— f language has allowed 
us to maintain a clean separation between interfaces and im- 
plementations (e.g., for the equation of state, Poisson equa- 
tion solver, and initial conditions etc.), and features such 
as templates have allowed us to write generic code without 
sacriflcing runtime performance. 

Since standard C-I--I- does not provide facilities for par- 
allel computing, it is necessary to use additional libraries to 
handle the parallelisation. We have chosen to implement par- 
allelisation via the Message Passing Interface (MPI). Since 
both optimising, ISO-compliant C++ compilers and high 
quality MPI implementations are available for virtually ev- 
ery major computing platform, our code is highly portable. 



4 SOLVING THE POISSON EQUATION 

Equation JSJ is distinct from equations 11131 in that it re- 
quires global, rather than local, information. There are a 
number of methods that can be used to solve the Pois- 
son equation. These include general elliptic equation set 
solvers, multigrid met hods, mult ipolc met hods, and spec- 
tral methods (see e g. iMotl et alJ.2002: Fryxefl et al.ll200d : 



iMuller fc Steinmet j[l995l: IStone fc Normanlll992^ . Spectral 
methods tend to be the most efficient, and implementing 
them on a regular Cartesian grid is straightforward. 

The solution of the Poisson equation requires the spec- 
ification of a boundary condition on some closed surface. In 
most physical problems, this surface is chosen to lie at infln- 
ity, upon which the potential is chosen to vanish. However, 
since our computational domain is flnite, it is not possible to 
impose a boundary condition at infinity in a straightforward 
manner. Instead, we define the value of the potential on the 
surface of our domain, which we compute via a multipole 
expansion: 



-r ( 



;=o m^-l 



21 + 1 



^Y'im(x) , 



where 



dx'r"Yil,{±')p{ji') . 



(25) 



(26) 



In practise, it is only necessary to include the first few multi- 
poles (for our purposes Zmax = 5) to obtain accurate bound- 
ary values. Note that the boundary condition at infinity is 
built into the multipole expansion. 

Given the Dirichlet boundary condition, it is possible to 
solve Poi sson equation via a discrete sine transform (DST) 
(see e.g. iPress et alJll992^ . Written in its finite-difference 
form, becomes 



E 



-I"-! 



2$ + $_ 



(Aa 



(27) 
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In terms of their discrete sine transforms "I> and p, $ 
and p are given by 

7-1 J-l K-l 



J. z \ - \ - \ - 5: . Tvim . ivjn . nkp 

" jjK 1^2^!^ ~r ~ ''''' 'w 

(28) 



m — 1 n — 1 p— 1 



7-1 J-l K-l 



UK 



Ev^ • Trim . ttto . Trfep 

2^ 2^ Pm.n.p sm — p sm — sm — 



m — l n—l p—l 



(29) 

where i, j, k, and /, J, K define the location in, and the 
dimensions of, the computational domain, respectively. Sub- 
stituting these expansions into 1271 gives 



-47rG 



Pm,n,p 



(30) 



where 



+ 



(Ax)2 

2 
2 



J 

cos-j 



(A^)2 

The potential is then computed from 12811 . 

Expanding "1> in terms of the sine basis functions of the 
Fourier series ensures that it vanishes at the boundaries of 
the domain. Non-zero boundary conditions can be incorpo- 
rated by adding an appropriate source term to the right side 
of equation 1271 1. We may define = "I> — "I>^ where now 
$^ is determined by equation 12511 at one zone beyond the 
boundary and vanishes everywhere else. The resulting equa- 
tion for is the same as equation 1271 in the interior and 
is given by 

E ^^-^S^-4.Gp-^^4.Gp',(31) 

i — x,y,z 



(Ax«)2 



on the ±jth boundary. As a result, the effective source terms 
are given by 

1 



(As)2 

1 

(A^ 
1 



(32) 



(A2)2 

To summarise, our procedure for solving the Poisson 
equation is: 

(i) Calculate <1>^ via the multipole expansion 1251 . 

(ii) Calculate the effective source terms for "!>' from 1321 . 

(iii) Perform a DST on the effective source terms. 

(iv) Calculate from CT . 

(v) Perform a DST on $' to determine 

We do not actually need to add $^ to our final answer since 
it only affects the ghost points outside our grid. Note that, 
because we use a second-order finite-difference to determine 
the gravitational acceleration in equation l|18^ . it is neces- 
sary to define $ on an extra surface of ghost cells on each 
edge of the domain. 
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Figure 2. A square pulse that has been advected five times its 
initial width (50 cells) using the donor cell (open circles), van Leer 
(filled triangles), and PPA (open squares) upwinding schemes. For 
reference, the original pulse profile is also shown. 



The DST is most efficiently parallelised in terms of a 
slab decomposition of the grid, as opposed to the ideal de- 
composition for the sourced advection step (which is cubi- 
cal). As a result, a significant amount of interprocess com- 
munication is required to prepare for the solution of the 
Poisson equation at each source sub-step. However, we have 
found that the time saved by using the DST more than out- 
weighs the penalty incurred by the communication overhead 
compared to alternative methods. 



5 TEST PROBLEMS 
5.1 Advection 

In order to test the advection scheme, we considered the 
advection of a square pulse (without source terms). In Fig- 
ure |51 the pulse is shown after being advected five times its 
initial width (50 cells) using both the donor cell and van 
Leer upwinding methods. It is clear that both methods are 
diffusive, with the donor cell method substantially more so. 

In general, diffusion will lead to errors in both the am- 
plitude and the phase of an advected pulse. In order to 
quantify these errors for diffusion resulting from the upwind- 
ing scheme, a sine wave was advected with periodic bound- 
ary conditions for 100 times its wavelength. By this time, 
the donor cell upwinding scheme has diffused the sine wave 
completely, hence only the van Leer and PPA methods are 
shown in Figure |21 The errors are at the 4% and 0.4% lev- 
els, respectively, with deviations becoming most significant 
at extrema. In both the square pulse and the sine wave, a 
noticeable asymmetry (which is determined by the direction 
of propagation) develops as a result of higher-order effects 
in the upwinding schemes. 
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Figure 3. A sine wave is advected with periodic boundary con- 
ditions for 100 times its wavelength (200 cells) using the van Leer 
(filled triangles) and PPA (open squares) upwinding schemes. In 
the top panel the density profile is shown explicitly, while in the 
bottom the residuals are plotted. For reference the analytical re- 
sult is also shown. 



5.2 Sod Shock Tube 

The pressure source term in equation Q was tested by the 
Sod shock tube problem. The Sod shock tube consists of an 
initial density and pressure discontinuity, and its subsequent 
evolution for an ideal gas (F = 1.4) and a specific set of 
initial conditions. For x > 0, p — 0.125 and P — 0.1, while 
for a; ^ 0, p = 1 and P = 1. Because it is the entropy density 
and not the pressure that is evolved, it is necessary to find 
s as a function of p and P for an ideal gas: 

' ' (33) 



s = In 



where n — 



,p"^^J F — 1 

The Sod shock tube is useful as a test because the result- 
ing p and P profiles for any given time can be ca lculated 
analytically (see, e.g.. Is^lioTa iHa.wlev et al.lll984^. 

In Figure |^ the numerical results from our code are 
compared to the analytical solutions. Overall, they are in 
good agreement, with the exception of two minor discrepan- 
cies. The most notable discrepancy is the entropy deficit in 
the post-shock fluid (0.184 < a; < 0.35). This is a result of 
using the adiabatic condition, and thus ignoring entropy pro- 
duction at shocks. Hence, the higher analytical value is easy 
to understand. Because we intend to apply the code to sce- 
narios in which the adiabatic condition holds to a very good 
approximation, we expect the entropy deficit to be physi- 
cally insignificant. The second discrepancy is the overshoots 
at points where the s lopes o f quantities change discontinu- 
ously. As discussed in lStone fc Norman (1992), this is a real 
result, originating from the numerical viscosity inherent in 
any finite-difference code. The most important result, how- 
ever, is the fact that the artificial viscosity causes the shock 
fronts to be well behaved in our code. 



Figure 4. The density, pressure, velocity, and entropy are shown 
for the Sod shock tube at i = 0.2 (the units of which depend upon 
the units chosen for the pressure and density). 200 cells were used 
with van Leer upwinding. The head and tail of the rarefaction 
wave are located a,t x = —0.235 and x = —0.014, respectively. The 
contact and shock discontinuities are at x = 0.184 and x = 0.35, 
respectively. 



5.3 Pressure- Free Collapse 

The gravitational source term in equation ^ was tested via 
the pressure-free collapse of a uniform density sphere. Once 
again, there is an analytical solution: 



r = ro cos (3 

-6 a 

p — po cos p 



/3+ -sin 2/3 



(34) 



-Gpo 



-1/2 



(see, e.g.. lStone fc NormanllOO^) . Figure|S]depicts the result 
after allowing the radius to halve (at t = 0.909 for G — 1), 
for a 256 x 256 x 256 cell grid. There is a small excess on 
the edges resulting from our implementations of viscosity 
and consistent transport (which necessarily treats the ad- 
vection of velocity into the edges differently due to the den- 
sity gradients). Overall, it does show good agreement with 
the analytical prediction. 



6 APPLICATION TO A PULSATING WHITE 
DWARF 

6.1 Hydrostatic Equilibrium 

The problem of choosing an equilibrium fluid configuration 
is made non-trivial by the finite differencing of the the dy- 
namical equations. Consequently, a method to produce an 
equilibrium solution for the finite difference equations is re- 
quired. For a barotropic equation of state, we have chosen to 
make use of the self-consistent field (SCF) method (see, e.g.. 
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Figure 5. The numerical (open circles) and analytical (solid 
line) solutions for the density as a function of distance along a 
radial section for the pressure free collapse of a uniform density 
sphere are shown. The initial radius and total mass of the sphere 
was unity. A 256 X 256 X 256 cell grid was used. With Newton's 
constant given by G = 1, this occurs at t = 0.909. 



iMotl et aljl2003 : lHachisulll986l : lOstriker fc Mark^ll968^ . Be- 
cause it is well described elsewhere, we will only summarise 
the procedure here. 

(i) An initial guess for the density (taken from the contin- 
uous solution) is used to generate the gravitational potential 
via the method described in 2] 

(ii) The new gravitational potential and the initial den- 
sity guess are then used to calculate the Bernoulli constant 
at the centre of the star. 

(iii) the Bernoulli constant and the new gravitational po- 
tential are used to calculate the enthalpy at all points on the 
grid, which is then subsequently inverted to yield the new 
density guess. 

This procedure is iterated until the Bernoulli constant con- 
verges to some specified tolerance-i.e. when the fractional 
change is less than some small value (say, lO"'^^). The re- 
sulting density distribution is a solution to 



h4 



h. 



$4 



and, hence, no net momentum flux is generated if the source 
terms are given by equation (1191 . Note that, if the source 
terms are given by equation l|18|l , this may still produce a net 
momentum flux, and is not necessarily a good approximation 
to equilibrium in that case. 
When 



dP 



P 



(36) 



Model 


M(Mq) 




uj, (Hz) 


ujf2 (Hz) 


ujp2 (Hz) 


CWD 


0.632 


8.56 


0.365 


0.562 


1.15 


HWD 


0.632 


11.2 


0.243 


0.560 


0.749 



Table 1. Stellar properties for a cold white dwarf with (CWD) 
and without (HWD) an isothermal envelope. Specifically, the 
mass, radius, fiducial stellar frequency oj, = i^/ GM/R^, frequency 
of the adiabatic quadrupolar fundamental mode, and the fre- 
quency of the lowest order adiabatic quadrupolar p-mode. Note 
that the inclusion of the isothermal envelope does not change the 
mass appreciably while significantly increasing the radius. 




the pressure gradient required to preserve hydrostatic equi- 
librium cannot be resolved on the grid. For a star, this can 



Figure 6. Shown in the top panel are the density profiles for the 
cold white dwarf with (solid) and without (dashed) the isother- 
mal envelope. The two lower panels are the radial displacement 
profiles for the quadrupolar fundamental mode (/2) and the low- 
est order quadrupolar p-mode (p2) for the two models. Note that 
the density and /2 mode profiles are very nearly the same for the 
two cases. However, the mode profiles differ substantially for the 
p2 mode. 



result in strong, inwardly directed forces at the surface, driv- 
ing shocks into the interior. We have found that adding an 
isothermal envelope can mitigate this problem by pushing 
the region where this inequality is true off the grid, while 
adding an insignificant amount of mass to the star itself. 
This is done explicitly by setting a fiducial density (which 
we chose to be 10~^ of the central density) at which the 
equation of state changes from that of a cold white dwarf 
to a r = 1 polytrope. The polytropic constant is chosen 
such that P{p) remains continuous across the transition. Ta- 
ble Q compares the properties of the cold white dwarf with 
(HWD) and without (CWD) the isothermal envelope. Note 
that while the isothermal envelope increases the radius sig- 
nificantly, it does not change the mass or the frequency of 
the quadrupolar fundamental mode {u}f2)- The reason for 
this can be seen in Figure il The /2 mode is more strongly 
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weighted in the core where the addition of the isothermal 
envelope makes no difference. In contrast, the lowest-order 
quadrupolar p-mode is substantially affected by the pres- 
ence of the envelope. This probably results from the fact 
that the radial wavelength of the p2 is much closer to the 
height of the isothermal envelope. Henceforth, all evolutions 
were begun with the HWD model listed in Tabled 

The quality of the equilibrium generated by the SCF 
method may be explicitly demonstrated. FigureQshows the 
evolution of the centre-of-mass position, net momentum, and 
the fraction of the total energy that is converted into kinetic 
energy for a star initially in hydrostatic equilibrium. The 
last quantity is given in terms of the kinetic, internal, and 
gravitational components: 

J^pv^d^x, n = Jp<fx, W ^ j p'^d^x. (37) 

Despite an initial exponential rise, these quantities saturate 
at relatively low levels for all resolutions shown. Note that 
all times are measured in dynamical times of the cold white 
dwarf, tcwD = which is approximately the time it 

takes for a disturbance to cross the star. 



6.2 Oscillation Modes 

In general, the problem of interest is dynamical. Specifically, 
we are interested in the non-linear evolution of the oscilla- 
tion modes of a cold white dwarf which are being excited 
resonantly by tidal forces. Towards this end, it is impor- 
tant to obtain a measure of the numerical quality factor (Q; 
the e-folding time of the energy in the oscillation), and the 
oscillation frequencies themselves. That the latter may be 
different from the frequencies in Table Q is a result of both 
the approximation of discrete cells and the fact that the 
finite-difference equations are distinct from the continuous 
equations. However, we expect the deviation to be small, and 
therefore a close agreement between the predicted and ob- 
served frequencies serves as yet another test for the correct- 
ness of our code. Both the quality factor and the oscillation 
eigenfrequencies can be obtained by deforming the star in a 
particular way, and analysing the subsequent oscillations. 

We deformed the star by adding a fractional quadrupo- 
lar perturbation to the density, i.e. 



/\p{v) = Ap{r)Y^^{6,d}), 



(38) 



where the amplitude. A, was chosen to be small (10~^) so 
that the resulting oscillation occurred in the linear regime. 
This initiated an even m — 2 standing wave on the star. 
Figures |H| and |^ show the resulting evolutions for a num- 
ber of grid resolutions. The same diagnostics as those used 
to demonstrate hydrostatic equilibrium are shown in Figure 
[7| In this case as well, the centre-of-mass and momentum 
drift saturate at levels well below those of interest. Unlike 
hydrostatic equilibrium, there now exists a non-vanishing 
kinetic energy. It is strongly harmonic and decays exponen- 
tially. Because the initial perturbation excited all of the even 
quadrupolar modes with m — 2, there are a number of dis- 
tinct decay constants, with the slowest being due to the /2 
mode. This exponential decay at late times may be fit to 
estimate the numerical Q, found here to be on the order of 
6000. 

In Figure 1^1 the quadrupolar moments are shown. The 
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Figure 7. Shown arc the centre-of-mass (top panels), net momen- 
tum (middle panels), and fraction of the total energy converted 
into kinetic energy (bottom panels) for a number of grid resolu- 
tions (note the different time scales). In all cases these quanti- 
ties saturate well below significant levels (e.g., for the worst case, 
the centre-of-mass moves by less than 10~* cell widths in the 
150 dynamical times shown, thus it would require roughly 10^^ 
dynamical times before the centre-of-mass moves one stellar ra- 
dius. Typically, these appear to turn over, implying that they may 
never rise significantly above 10~^ cell widths.) 
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Figure 8. Same as Figure |7| for the case when a quadrupolar 
perturbation is present (note the difference in scales in comparison 
to that figure). The white line drawn through the oscillations is 
for a Q of approximately 6000. 



Figure 9. The quadrupolar moments of the perturbed star for 
each of the resolutions considered in Figure |S] From top to bot- 
tom, the panels are the odd m = 2, odd m = 1, m = 0, even 
m = 1, and even m = 2 moments. Note the difference in scales 
of the different moments, namely that the even m = 2 moment is 
two orders of magnitude larger than the m = moment (which 
is sourced by the Cartesian geometry) and nine orders of magni- 
tudes larger than the others. 
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Figure 10. Shown are the power spectra of the even m = 2 
quadrupolar moment as a function of angular frequency (using the 
mean squared amplitude normahsation). As expected, for each 
grid resolution there is a strong spike coincident with the /2 mode 
frequency predicted for the HWD model. 



even m = 2 moment is strongly dominant as expected. 
It also has a very clear harmonic structure. This may be 
Fourier analysed to produce the dominant oscillation mode, 
as shown in Figure [TUl In the power spectrum of the even 
m = 2 quadrupolar moment, there is a peak which extends 
five orders of magnitude above the rest of the spectrum. This 
peak is clearly identifiable with the /2 mode, and appears 
to have very nearly the frequency predicted by the HWD 
model. 



7 CONCLUSIONS 

We have developed and tested a parallel, simple and fast 
hydrodynamics code for multi-dimensional, self-gravitating, 
adiabatic flows. Both the advection terms and the solution 
for the self-gravity are greatly simplified by the use of a 
uniform Cartesian geometry, ultimately leading to explicit 
conservation of mass, entropy, and momentum to nearly nu- 
merical accuracy. The simplifying assumption of adiabaticity 
and the absence of shocks eliminate the necessity for more 
numerically expensive schemes, yielding an efficient code. 

We have also applied our code to a number of standard 
diagnostic problems in order to verify its physical correctness 
and limitations. This has been done in a systematic fashion 
intended to test each aspect of the code separately, includ- 
ing the advection scheme, the pressure source terms, and 
finally the gravitational potential. Finally we have demon- 
strated the fitness of the code for the problem which mo- 
tivated its development: the study of tidally excited adi- 
abatic oscillations on white dwarfs. This has been done in 
two stages: firstly, verifying the long term numerical stability 
of a white dwarf in hydrostatic equilibrium, and, secondly. 



measuring the numerical quality factor (found to be on the 
order of 6000) and the quadrupolar fundamental mode fre- 
quency (found to be very nearly that predicted by a linear 
mode analysis of the white dwarf model). The details of 
tidally exciting adiabatic oscillations with non-linear ampli- 
tudes on cold white dwarfs, and their subsequent evolution 
will be discussed in a future publication. 
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